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Abstract 

In the light of recent experimental findings that gap junctions are 
essential for low level intensity detection in the sensory periphery, the 
Greenberg-Hastings cellular automaton is employed to model the re- 
sponse of a two-dimensional sensory network to external stimuli. We 
show that excitable elements (sensory neurons) that have a small dy- 
namical range are shown to give rise to a collective large dynamical 
range. Therefore the network transfer (gain) function (which is Hill 
or Stevens law-like) is an emergent property generated from a pool 
of small dynamical range cells, providing a basis for a "neural psy- 
chophysics" . The growth of the dynamical range with the system size 
is approximately logarithmic, suggesting a functional role for electrical 
coupling. For a fixed number of neurons, the dynamical range displays 
a maximum as a function of the refractory period, which suggests ex- 
perimental tests for the model. A biological application to ephaptic 
interactions in olfactory nerve fascicles is proposed. 
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1 Introduction 



Experimental results show that individual receptor neurons have a 
small dynamical range ^ |2] • Given a properly defined stimulus inten- 
sity r (e.g. light intensity for rods in the retina, or odorant concen- 
tration for sensory neurons in the olfactory epithelium), the response 
f(r) of a single neuron (e.g. the spike frequency) saturates at a value 
fmax for sufficiently large values of r. Defining no and rgo as the 
stimulus values for which the response attain 10% and 90% of fmax, 
the dynamical range Ar is usually defined as 



in log units. In other words, the dynamical range roughly measures 
the number of decades over which the stimulus can be properly dis- 
criminated 1 , discarding stimuli too faint to be considered "detected" 
(r < rio) or too close to saturation (r > rgo). In olfactory receptor 
cells of the tiger salamander |T and the frog [2], Ar ~ 1 log unit 



The psychophysics literature, on the other hand, shows that our 
ability to discriminate external signals of different kinds covers a large 
dynamical range [H] ■ This is reflected in phenomenological psychophys- 
ical laws used to fit data, which state that the psychological percep- 
tion of a given stimulus r increases as ~ r a , a < 1 (Stevens' law), 
as ~ log(r) (Weber-Fechner law), or as ~ r a /(c a + r a ) (Hill func- 
tion, which may be thought of as a saturating Stevens' law). That 
is, differently from individual sensory neurons, organisms perform sig- 
nal compression, being able to process environmental stimuli whose 
intensities usually span several orders of magnitude. 

How can we reconcile these results? How could this signal com- 
pression be physically implemented considering that, at the very first 
stage of the signal processing, individual sensory neurons have a small 
dynamical range? 

Following a conjecture put forward by Stevens |3], we propose that 
this signal compression occurs already at the periphery of the nervous 
system, that is, at the first "relevant" synaptic level. In our model, 

1 An alternative definition employs 5% and 95% of the saturation value for evaluating 
Ar. This is of course just a matter of choice. Our results remain qualitatively the same 
in either case. 





(10 dB). 
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this processing comes out naturally as a result of the collective phe- 
nomenon of non-linear amplification in excitable networks. While 
individual sensory neurons have a small dynamical range, an array 
of connected neurons can collectively exhibit a very large dynamical 
range. 

The motivation for the model comes from recent experimental find- 
ings showing that electrical synapses by gap junctions are present 
in the periphery of different sensory systems, such as the olfactory 
glomeruli [3J [3] and the retina [Bj. With this information at hand, 
we model a given sensory periphery as a two-dimensional array of ex- 
citable elements which are laterally connected by electrical synapses. 
Experimentally, the topology of the neural network connected by gap 
junctions is still not known, and specific details probably depend 
on the particular sensory system under consideration. So the two- 
dimensional lattice can be considered as the natural first step towards 
a biologically acceptable model. As a possible biological application, 
it could represent the cross section of a nerve fascicle containing highly 
packed unmyelinated axons coupled by electric interactions, which tar- 
get a common secondary neuron, as found in the olfactory periphery. 
The one-dimensional case, which may be thought of as an approxi- 
mation for networks where neurons are coupled to two neighbors in 
average, has been previously addressed in Ref. |7j. 

In section |2] we introduce the cellular automaton model used for 
each excitable element and briefly review results from the literature. 
Our results are discussed in section while section 0] brings our con- 
cluding remarks and suggestions for experimental tests. 

2 The model 

In order to simulate large networks, we make use of the simplest pos- 
sible model for an excitable element: an n-state cellular automaton 
(CA). The state of the i-th. cell (i = 1, . . . , N) at discrete time t is 
denoted by Xi(t) 6 {0, 1, . . . , n — 1}. A field hi(t) £ {0, 1} indicates 
whether the stimulus at site i is supra- or sub-threshold. The CA 
model contains two ingredients: 

1. a cell spikes in the next time step (xj(t + 1) = 1) only if it is 
currently stimulated (hi(t) = 1) AND in its resting state (xi(t) = 
0) 

2. after a spike, a refractory period takes place (xj = 2, 3, . . . , n— 1), 
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during which no further spikes occur, until the cell returns to its 
resting state (formally, if Xi(t) > 0, then Xi(t + 1) = [xi(t) + 1] 
mod n with probability one) 

This corresponds to one of the several variants of the model pro- 
posed by Greenberg and Hastings [Hj. It was initially designed to 
provide a simple explanation for the mechanisms underlying pattern 
formation in excitable media. Indeed, the onset of self-sustained or 
re-entrant activity (with e.g. spiral waves) in a plethora of systems (in- 
cluding from heart tissue to chemical reactions) has been a research 
topic of considerable interest (see e.g. Refs. and references 
therein). Different flavors of the Greenb erg-Hastings cellular automa- 
ton (GHCA) have been used in a variety of applications: collective os- 
cillations of pyramidal cells in the hippocampus |11U12| : noise- induced 
memory ^3] and the evolution of the HIV infection ^1]. In this paper 
we address an apparently unaccounted feature of the widely studied 
GHCA: the collective enhancement of the dynamical range of the sin- 
gle cells. 

We model the response of a sensory neural network to some exter- 
nal input, and assume that this response vanishes if input is removed. 
This is accomplished with the rules defined above which, together with 
a resting initial condition (xi(t = 0) = 0, Vi) and the synchronous up- 
date of the lattice, prevent self-sustained activity from occurring |12j . 
We also assume that the input signal h{t) (which may have been pre- 
processed by earlier layers in a biological network) arriving on cell 
i at time t is a Poisson process of supra-threshold events of stereo- 



typed unit amplitude. That is, Ii(t) = Yl n ^ (Mri ) where 5(a,b) is 



the Kronecker delta and the time intervals t n+l — tn are distributed 
exponentially with average (input rate) r, measured in events per sec- 
ond. At each time step, the probability of arrival (per neuron) of an 
external input is therefore 



where r = 1 ms coincides with the approximate duration of a spike 
and is the time scale adopted for a time step of the model. With this 
choice, the number of states n of the GHCA corresponds roughly to 
the refractory period (measured in ms) of the sensory neurons (or the 
sensory axons). 

The network connectivity is implicitly defined by the field hi. For 
uncoupled cells, we have simply h^\t) = Ii(t), i.e. cells are excited 




A( r ) = 1 - e 



—rr 



(2) 
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only by external stimuli. We model the coupling via gap junctions 
with the field 



hf\t) = l-[l-I i (t)]H[l-5(x,(t),l)] 



(3) 



j 



i.e. cells are excited by external stimuli or at least one spiking cell j 
(where j runs over the neighborhood of i). We have employed periodic 
as well as open boundary conditions, with indistinguishable results. 
The number of sites in the regular lattice with coordination number 
z is N = N x N y . 

The response F(r) of the lattice to a given stimulus intensity r 
is given by the mean firing rate per cell in a long time interval T 
(typically we have used T = 0(1O 5 ) time steps). This should be 
compared to the mean firing rate per cell / of an ensemble of uncoupled 
cells, whose dependence on r can be exactly calculated: 



This linear saturating function will be the benchmark against which 
results from the lattices should be compared. Note that, with the 
definition of Eq. ^ uncoupled cells have a dynamical range 



which is a smooth curve and, for moderate values of n, is well ap- 
proximated by the asymptotic value Ar ~ log 10 (81) ~ 19 dB. This 
dynamical range is slightly larger than the experimental results, but 
still too small to account for psychophysical data [3]. 



For very low values of r, the stimulation of a single site generates a 
wave of activation which propagates through the entire lattice (see 
Fig. ^a)). Before another stimulus is produced, the wave front is 
annihilated at the walls (with open boundary conditions) or when it 
touches itself (with periodic boundary conditions), because sites which 
have spiked in the last n — 1 time steps cannot spike again due to the 







3 Results and Discussion 
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Figure 1: Snapshots of a triangular 150 x 50 lattice where 10— state CA 
interact with their z = 6 nearest neighbors (periodic boundary conditions). 
Closed circles: Xi(i) = 1 (spiking). Open circles: 2 < Xi(t) < n (refractory). 
Background: Xi(t) = (resting). Panel (a): r = 0.005 s _1 . Panel (b): 



refractory period. In this regime, a single stimulus produces N spikes, 
yielding an TV-fold amplification, 

F r ~° NX ~ Nrr . (6) 

At the other extreme, for very high values of r the propagation 
of waves is suppressed by excess of external stimulus: in this regime, 
the coupling becomes irrelevant since each cell is already spiking at 
its maximal rate, 

p r ~°° t ~ _ = f (7) 

A — J — — - 1 max ) (' J 

n 

yielding no amplification, F/f ~ 1. 

In between those two regimes, waves are created, propagate and 
annihilate when meeting one another, yielding complex geometric pat- 
terns which appear from the remains of partially destroyed waves (see 
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Fig. IHb)). The interplay between wave creation and wave annihila- 
tion provides the self-limited mechanism by which the amplification 
F/f decreases smoothly from O(N) (for small r) to 0(1) (for large 
r). This is a mechanism which leads to signal compression. 

Fig. Efa) shows the result of numerical simulations for square lat- 
tices with z = 8 and n = 10 in a log-log plot. The linear response 
(Eq. EI) for low values of r can be clearly seen at the leftmost portion 
of the graph, requiring lower values of r for larger system sizes (for 
a fixed low value of r, the larger the system the more probable it is 
for two excitable waves to interact). As r increases, wave annihilation 
starts competing with wave creation, bending down the nonlinear re- 
sponse. Finally, saturation ensues for values of r close to 1/n, making 
the coupled and uncoupled cases nearly indistinguishable. The en- 
hancement of the dynamical range becomes clear when we plot F(r) 
in a linear-log plot [Fig. Efb)]: since the amplification due to coupling 
is larger for lower r, the curves stretch upward, compressing more 
decibels of stimuli in the same [O.l-F^axi 0*9-^max] response interval. 

The dimensionality of the lattice has a significant impact on the 
collective response, as exemplified by the results for a one-dimensional 
chain in Fig. |5J The differences can be understood as consequences 
of the following factors: a) in ID wave fronts have a constant size (2 
sites), while in 2D their size grows with time; b) because of a), waves 
in 2D can be partially destroyed (as exemplified in Fig. 1), while 
in ID they are completely annihilated upon collision; c) the time it 
takes for a wave to reach the borders scales with N in ID, but with 
y/~N in 2D. The combination of factors (a), (b) and (c) gives rise to 
finite size effects which are much stronger in 2D than in ID. As Fig. 13 
shows, for example, ID chains with 1600 or 400 neurons would have 
essentially the same dynamical range, since finite size effects appear 
below rio = 10 Hz. For 2D lattices, on the other hand, the dynamical 
range changes. 

We have simulated two-dimensional square lattices with z = 4 and 
z = 8 neighbors, as well as triangular lattices with z = 6. For a 
fixed dimensionality, the differences in coordination number are not 
important to our investigations, as can be seen in Fig. |2] (the curve 
for z = 6 falls between those of z = 4 and z = 8). In what follows, we 
focus on results for z = 8. 

For a given refractory period n, the size of the lattice regulates the 
crossover from the linear to the nonlinear response. If the sensitivity 
level rio lies within the linear regime (as it does, for instance, in Fig.|2J), 
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Figure 2: Response curves for square lattices of 10— state automata, (a) 
coupled (F, filled symbols) and uncoupled (/, open symbols) cases for z = 8 
neighbors. The dashed line corresponds to the exact solution (Eq. while 
the dotted lines are the linear approximation, Eq. |U1 Open triangles are 
results for a one-dimensional coupled system; (b) Linear-log version of (a), 
with the dot-dashed line representing a square lattice with N = 40 x 40 and 
z = 4 neighbors. Arrows indicate the dynamical range. 



then the dynamical range will be affected by the system size. It is 
important to note that this is not an artificial feature due to the 
somewhat arbitrary definition of Ar. Any reasonable definition of 
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Figure 3: Response curves of 3— state (filled diamonds) and 10— state (open 
diamonds) automata for different lattice sizes iV (N = 25 for n = 3 and 10, 
up to N = 51529 for n = 3). The upper (lower) horizontal line corresponds 
to 90% (10%) of the saturation firing rate for n = 3 (F max = 1/3 ms _1 ~ 
333 s _1 , thick dashed) and n = 10 (F max = 1/10 ms _1 ~ 100 s -1 , thin 
dotted). Inset: sensitivity (r 10 , circles) and saturation (r 90 , squares) levels 
(filled symbols for n = 3, open symbols for n = 10) as a function of the 
number of sites N in the lattice. The vertical arrow illustrates the dynamical 
range for a 40 x 40 lattice with n = 3. Dashed lines correspond to Eq. |H1 



the dynamical range will reveal a sensitivity gain (decrease in rio) for 
increasing values of N, as depicted in the family of curves of Fig. |31 
Note that the saturation level rgo also changes, but much less than 
rio, leading to an enhancement of the dynamical range. 

It is also worth mentioning that, since rio an d rgo are defined only 
relative to a saturation response, these quantities will depend on n. 
Fig. |31 and its inset illustrate this fact for n = 3 (filled symbols) and 
n = 10 (open symbols). Note that the relative meaning of "sensitivity" 
and "saturation" explains the following apparent contradiction: even 
though the n = 3 curves lie above the n = 10 curves (for the same 
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system size), their sensitivity and saturation levels also lie above those 
of the n = 10 case. In other words, if rio is assumed to lie in the linear 
regime, then Eqs. Inland H3 yield a decreasing dependence on iV and n: 

1 / 0.1 \ 
r 10 ~--ln 1-— . (8) 



r V nN , 

This approximation is acceptable for moderate values of N (above 
which the interference among excitable waves invalidates the linearity 
assumption) and improves for increasing n (since this lowers r±o into 
the linear region — see inset of Fig. |2J). 

The above discussion only emphasizes that the dynamical range 
stands out as a better quantity to look at than no and rgo, being a 
bona fide dimensionless variable which could be compared to exper- 
imental data. The dependence of Ar on the system size can be seen 
in Fig. 0] Note that an expansion in Eq. |H]for N S> 0.1/n yields 

Ar(iV) ~ log 10 (r 90 ) + log 10 (10miV) . (9) 
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If we neglect the (weak) dependence of rgo on N, this expression ex- 
plains the approximately logarithmic growth in Fig. 0]before the curves 
level off. For larger values of N, the linearity assumption for rio breaks 
down and Eqs. |H]and|5]are no longer valid. 

In Fig. |I]one observes that for low values of N (say, up to 0(1O 3 )), 
lower values of n imply larger dynamical ranges (this cannot be in- 
ferred from Eq. since it omits the unknown and strong dependence 
of rgo on n). As N grows, however, the curves in Fig. 0] start leveling 
off for increasing values of n. This simple fact has very interesting 
consequences. Suppose that we fix the system size at one of the four 
values specified by the vertical lines in Fig. If we now vary the re- 
fractory period n, we obtain the four curves displayed in Fig. |S1 which 
have a maximum. 

If we accept the reasoning that a large dynamical range is a de- 
sirable property for a live organism, then the results presented so far 
suggest two biological outcomes: first, that it would be advantageous 
for the organism to have as many sensory cells as possible. By elec- 
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trically coupling them, an increase of more than 100% in dynamical 
range could be obtained, as compared to isolated neurons (see in- 
set of Fig. [5J) . This is supported by the experimental evidence that 
gap junctions are present in the sensory periphery (which was in fact 
what motivated the model 7 ) . Recent experiments in the mammalian 
retina [0] show that knocking out the gene responsible for connexin-36 
(which accounts for the neuron- neuron gap junction channels) leads to 
a dramatic decrease in the dynamical range of ganglion cell responses. 

However, the number of receptor neurons in a given sensory sub- 
system is usually limited by several factors, from metabolic costs to 
sheer occupied space (for instance, each olfactory glomerulus of the 
rabbit receives input from O(10 4 ) sensory neurons JH]). This leads to 
the second biological suggestion. Fig. |S] suggests that for a given fixed 
size of the subsystem, it would be advantageous for an organism to set 
the refractory period of the receptor neurons to the value n max that 
yields the maximum dynamical range. Since n max (N) seems to be 
an increasing function, we could expect the size of connected neuron 
patches and the refractory period of the neurons (which classically 
are two independent parameters) to be correlated. In other words, 
smaller patches would tend to have neurons with shorter refractory 
period. Naturally, it is important to check the robustness of these 
results in other models of excitable media before specific experimental 
tests could be proposed. 

4 Concluding Remarks 

We have presented results concerning the nonlinear response of the 
two-dimensional GHCA to a Poisson stimulus. In particular, the 
coupling between the excitable elements has been shown to greatly 
enhance the dynamical range due to the collective phenomenon of 
self- limited amplification of excitable waves. 

We stress that we do not claim this is the only explanation for 
the enhancement of the dynamical range. Variation of the values of 
rio and rgo within the population of receptor neurons, for instance, is 
another hypothesis. In that scenario, different stimulus ranges would 
"recruit" different groups of neurons, therefore leading to an overall 
enhancement of the dynamical range (even if the neurons were uncou- 
pled). Theoretical work in this direction by Cleland and Linster |16j 
has shown that this could in principle be achieved by heterogeneous 
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overexpression of odorant receptors in the olfactory epithelium. Note, 
however, that doubling the dynamical range of the ensemble (as com- 
pared to that of individual neurons) would require the order of a 
hundred-fold overexpression of receptors, in that model. As Cleland 
and Linster point out, "measured spare receptor capacities in intra- 
cortical and culture systems studied to date are typically less than 
twofold" ng. 

Interestingly, the different proposed mechanisms are not mutually 
exclusive. In fact, "recruitment" as described above could in princi- 
ple coexist with the self-limited amplification we propose. And both 
could cooperate with further nonlinear post-processing of the signals 
(e.g. by mitral cell integration in the olfactory bulb). We believe that 
the strength of our model lies on two central issues. The first one is 
experimental evidence: the model suggests a functional role for the 
electrical synapses by gap junctions which have been experimentally 
found in the sensory periphery. The second point is that the mech- 
anism is simple. Indeed, the model relies more generally on lateral 
excitation, and could in principle even forego gap junctions specifi- 
cally. In the olfactory nerve, for instance, it could be implemented via 
the electrical coupling mediated by ephaptic interactions among neigh- 
boring axons, which have recently been modeled by Bokil et al. |17j . 
Neuronal circuits where mutual bidirectional excitation by chemical 
synapses occur could also be subjected to a similar modeling. 

The general question of how brains represent sensory stimuli is a 
classical and very difficult one, and has been addressed by physicists 
like Fechner, Plateau and Maxwell. Our model addresses only part of 
the problem, containing a mechanism to enhance the coding of a very 
basic (perhaps the most basic) property of the senses: intensity of the 
stimulus. 

Consider, for instance, the retina. Our primary concern is with 
electrical coupling not at networks of non-spiking neurons (as cones, 
horizontal cells etc) but at excitable networks of spiking cells. So 
we are interested in wave generation and synchronization at the first 
relevant spiking network (ganglion cell layer, in the case of vision) 
as a candidate mechanism for intensity coding in the output of optic 
nerve fibers. So, the "input" level r presented in our model should be 
interpreted as the input to the appropriate excitable network which 
we speculate to be situated at the ganglion level, not the raw in- 
put to the retina. There is experimental evidence that the spikes of 
ganglion cells are correlated with small latency (typical of electrical 
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coupling) |181 119j . which is consistent with our model. For specificity 
coding, on the other hand, other mechanisms seem best suited, such as 
lateral inhibition (which is well established in the literature, occurring 
in previous layers |2U|). The olfactory epithelium, on the other hand, 
has no known topographical structure in the expression of different re- 
ceptors, presenting "mosaic" receptor expression In this case, the 
excitable elements of the model could be thought of as unmielynated 
axons inside an olfactory nerve fascicle (which are known to pertain 
to a single class of receptors and innervate a single glomerulus). Thus 
the model gives a potential explanation to the large dynamical range 
observed in the intensity coding of individual glomeruli, but not to 
the much more difficult problem of how the identity of different smells 
are combinatorially coded by the several glomeruli. 

Real neurons are obviously much more complex than the our simple 
CA model: refractory periods are not necessarily absolute, responses 
can come in form of bursts, adaptive firing etc. However, we have 
also examined more biologically realistic models previously, leading to 
results similar to those of the GHCA [7ll22|. The robustness of the re- 
sults presented in this contribution with respect to different excitable 
models, lattice connectivity and stimulus statistics is the subject of 
ongoing research. For instance, non-homogeneous two-dimensional 
lattices have been known to produce self-sustained activity (spiral 
waves) since the seminal work of Wiener and Rosenblueth 23 j. It 
remains to be investigated how to connect excitable elements in order 
to maximize the dynamical range of the collective response, but also 
suppressing spiral waves. 

Despite its lack of biological realism, the CA approach is extremely 
useful since it allows the simulation of large networks and even the 
calculation of collective properties under appropriate approximations. 
Indeed, analytical tools such as hierarchical mean field approxima- 
tions |24| are available for studying cellular automata. We have ob- 
tained good results for the one-dimensional system with a mean field 
approximation at the pair level [23], which analytically reveals the 
low-r signal compression. The single site mean field approximation 
for the GHCA fails in hypercubic lattices of any dimension, and we 
are currently investigating the next hierarchical step to analyze the 
two-dimensional system. It is interesting to point out, however, that 
unless the mean field calculation contains finite size corrections, it 
may not be the best tool for understanding the phenomena occurring 
in finite cell assemblies. As has been previously noticed regarding the 
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olfactory nerve fascicles, many interesting biological phenomena occur 
precisely at intermediate values of the lattice size, which emphasizes 
the usefulness of computer simulations in addiction to analytical re- 
sults that must assume large N (thermodinamic) limits. 
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